HEAD
El paso de alineación o mapeo de secuencias cortas reads producto de la secuenciación NGS a un genoma de referencia se considera una parte integral de los análisis de genomas (Pham-Quoc, Kieu-Do, & Thinh, 2019).
El problema del mapeo de reads consiste en ubicar o alinear dentro de un genoma de referencia previamente secuenciado los millones de reads para luego ensamblarlas a un nuevo genoma (Bak, Bodziony, Migdałek, Pareek, & Żukowski§, 2020) o en el caso de que no exista un genoma de referencia realizar un ensamble de novo.
Esta carga de trabajo se la realiza con diferentes enfoques de alineación, sin embargo, al igual que cualquier software, cada enfoque presenta diferentes formas de ejecución y consecuentemente de análisis.
Para realizar esta tarea se han desarrollado varios programas como lo son, BWA (Li & Durbin, 2009), Bowtie 2 (JHU, 2016) y otros, cuya función principal es hacer coincidir secuencias cortas provenientes de NGS con secuencias de referencia (Bak, Bodziony, Migdałek, Pareek, & Żukowski§, 2020). Actualmente, BWA es un software muy popular para fines de alineamiento, debido a su compatibilidad con datos de Illumina y la integración de distintos algoritmos de alineamiento con mayor rendimiento para reads de distintas longitudes (Li & Durbin, 2009).
El resultado del mapeo o alineamiento se guarda en un archivo d etexto plano con formato SAM (Sequence Alignment Map) o su equivalente comprimido en binario (BAM). Este archivo contiene de forma codificada toda la información de dónde y cómo cada read se alineo al genoma referencia. De el se pueden luego extraer información de variantes para análisis poblacional o estudios de asociación genómica o información de conteo de reads en regiones codificantes para en análisis de expresión de genes.
El formato SAM consta de un encabezado que comienza con el símbolo @ y una sección de alineamiento que contiene la información de cada uno de los reads que alineo al genoma de referencia. Para conocer más acerca del formato SAM puedes revisar el paper que describe el formato o simplmente en wikipedia. Los archivos SAM se pueden analizar y editar con el software SAMtools. A continuación se muestra un ejemplo del formato SAM con los 11 campos obligatorios del alineamiento y 1 campo opcional.
Formato SAM
BWA
Es un algoritmo de alineación para alinear lecturas de secuencia o secuencias de consulta largas con un genoma de referencia grande como el humano. Elige automáticamente entre alineaciones locales y de extremo a extremo y admite lecturas de extremos emparejados (Li, 2013).
Puedes visualizar las opciones y más información de bwa Aquí
Samtools
Es un conjunto de utilidades que manipulan alineaciones en los formatos SAM (Sequence Alignment Map), BAM y CRAM. Convierte entre los formatos, clasifica, fusiona e indexa, y puede recuperar lecturas en cualquier región rápidamente, puedes tener más información sobre este programa Aquí
Para facilitar el proceso solo se utilizará una muestra la que provieme de la base de datos SRA del NCBI y corresponden a lecturas crudas del salmón del Atlántico Salmo salar en formato fastq, obtenidas por secuenciación de extremos emparejados con un secuenciador Illumina HiSeq2000.
La muestra se alineará al genoma de la mitocondria del salmón del Atlántico puesto que dado que su pequeño tamaño, el alineamiento no debería tardar más de unos pocos minutos.
Descargar Genoma de referencia de la mitocondria e indexar.
Alinear muestra con genoma de referencia.
Analizar el resultado del alineamiento.
Conecte al servidor pomeo usando la siguiente dirección IP: 200.54.220.141 puerto 22.
Para Windows puede conectarse con el software PuTTy y las claves de acceso dadas en la primera clase.
Para configurar el canal bioconda se debe ejecutar el siguiente comando
conda config --add channels bioconda
Instalar los software bwa
conda install -c bioconda bwa
Instalar samtools requiere un proceso de instalación un poco más largo, siga cada paso con atención.
conda install -c bioconda samtools
conda config --add channels bioconda
conda config --add channels conda-forge
conda install samtools==1.11
Si desea constatar la ruta de instalación de un programa determinado puedes colocar el comando “whereis” en la terminal más el programa del que desea obtener la ruta.
whereis sratoolkit
whereis samtools
whereis bwa
La salida de cada comando indicará la ruta de instalación como sigue, para cada programa del que desees obtener información.
/home2/usuario/miniconda3/bin/bwa
Considerando que nuestras secuencias fastq originales tuvieron muy buena calidad, puedes trabajar directamente con ellas. Para ambos casos debes trasladar los archivos a una nueva carpeta, así:
Crear una carpeta denominada en tu usuario de home2
mkdir alineamiento
Ingresa a tu carpeta “alineamiento” y transfiere los archivos de clase anterior colocando los siguientes comandos en la terminal:
mv /home2/usuario/SRA_samples/SRR2006763/SRR2006763_1.fastq /home2/usuario/alineamiento/
mv /home2/usuario/SRA_samples/SRR2006763/SRR2006763_2.fastq /home2/usuario/alineamiento/
Lista tu carpeta de alineamiento para verificar que tienes lo necesario para el alineamiento, hasta ahora deben constar tus dos secuencias “SRR2006763_1.fastq” y “SRR2006763_2.fastq”
ls
Ahora descargaremos el genoma de referencia de la mitocondria de Salmo salar en la misma carpeta de alineamiento, para ello:
En tu navegador entra al siguiente link
Encontrarás una tabla con el listado del genoma de referencia donde se incluyen todos los cromosomas y el genoma de la mitocondria y a cuyo identificardor RefSeq darás clic, como se ve en la siguiente imagen
Listado genoma Salmo salar
Dar clic en la opción FASTA localizada bajo el título e identificador RefSeq de la referencia
Fasta option
Enviar la secuencia FASTA del genoma mitocondrial a un archivo
Genoma mitocondrial
Ve al lugar de descargas de tu computador y encontrarás el archivo denominado “sequence.fasta” al que debes renombrar como “mt.fasta” y subirlo a POMEO, para ello puedes trabajar con WINSCP
Si trabajas con WINSCP debes iniciar sesión, guardar tus datos y conectarte como se indica en la siguiente imagen
WINSCP
Al conectar debes ingresar nuevamente tu contraseña y aparecerá la interfaz de tu servidor con las carpetas que tienes creadas, aquí ingresarás en tu carpeta de alineamiento y arrastrarás el archivo descagado del genoma mitocondrial a la misma.
Alternativamente puedes cargar el genoma iniciando sesión online en el servidor POMEO a traves del puerto 8787 precargado con Rstudio server.
http://200.54.220.141:8787/auth-sign-in
Cuando termines lo anterior puedes ingresar a POMEO, listar tu carpeta de alineamiento con ls y verás tu archivo “mt.fasta” junto con tus secuencias fastq.
Una vez que incluiste en tu carpeta de alineamiento todos los archivos descritos en pasos anteriores podemos proceder a la etapa inicial del alineamiento, que corresponde a la indexación del genoma de referencia con bwa usando el siguiente comando:
bwa index mt.fasta
La salida del comando dará como resultado 5 archivos con extensiones “amb”,“ann”,“bwt”,“pac” y “sa”
Para el alineamiento tendremos las siguientes etapas
Alineamiento de las secuencias contra el genoma de referencia, cuya salida será un archivo con extensión “.sam”
Conversión del archivo SAM a BAM
Inspeccionar el archivo .sam de salida
Ordenar lecturas alineadas por posición
Indexación con Samtools
Exploración de datos con Samtools
Para ejecutar todas las etapas anteriores en ese orden se debe crear un script con nano denominado aln_mt.sh
nano aln_mt.sh
En el script coloca las siguientes instrucciones
#!/bin/bash -l
# para alinear tus dos secuencias fastq al genoma mitocondrial
bwa mem mt.fasta SRR2006763_1.fastq SRR2006763_2.fastq > SRR2006763.sam
# Transformar tu archivo sam a bam
samtools view -Sb -q 30 SRR2006763.sam > SRR2006763.bam
# ordenar tu archivo binario bam
samtools sort SRR2006763.bam -o SRR2006763.sort.bam
# indexar tu archivo bam
samtools index SRR2006763.sort.bam
Ejecuta tu script con bash
bash aln_mt.sh
Ahora que ya tienes tus archivos SAM/BAM puedes observar tu archivo sam con el comando less de linux (recuerda que es un archivo de texto plano)
less SRR2006763.sam
también puedes realizar un análisis estadístico estandar con los siguientes comandos
samtools flagstat SRR2006763.bam > muestra_stat.txt
Investiga otros comandos de samtools para explorar el alineamiento.
Algunas observaciones adicionales
BWA-MEM se recomienda para secuencias de más de ~70 pb, ya que, generalmente BWA-MEM es más tolerante con los errores dadas las secuencias más largas.
Si quieres ver lo que sucede en cada etapa de alineamiento, puedes colocar en la terminal cada comando del script de forma individual, e.g.
bwa mem mt.fasta SRR2006763_1.fastq SRR2006763_2.fastq > SRR2006763.sam # para alinear tus dos secuencias fastq al genoma mitocondrial samtools view -Sb -q 30 SRR2006763.sam > SRR2006763.bam # Transformar tu archivo sam a bam samtools sort SRR2006763.bam -o SRR2006763.sort.bam # ordenar tu archivo binario bam samtools index SRR2006763.sort.bam # indexar tu archivo bam
El archivo bam no se puede visualizar directamente porque está comprimido y en binario.
Procura listar tu carpeta cada vez que realices un trabajo, una descarga, o ejecutes un script para comprobar tus salidas.
Si quieres comprobar los tamaños de tus archivos al mismo tiempo que los listas, puedes hacerlo con
ls -l -h
Pham-Quoc, C., Kieu-Do, B., & Thinh, T. N. (2019). A high-performance FPGA-based BWA-MEM DNA sequence alignment. Wiley, 1-12.
Li, H., & Durbin, R. (2009). Fast and accurate short read alignment with Burrows-Wheeler Transform. Bioinformatics, 1754-60.
JHU. (2016). Johns Hopkins University. Obtenido de Manual of Bowtie 2 - Fast and sensitive read alignment: http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml
El paso de alineación o mapeo de secuencias cortas reads producto de la secuenciación NGS a un genoma de referencia se considera una parte integral de los análisis de genomas (Pham-Quoc, Kieu-Do, & Thinh, 2019).
El problema del mapeo de reads consiste en ubicar o alinear dentro de un genoma de referencia previamente secuenciado los millones de reads para luego ensamblarlas a un nuevo genoma (Bak, Bodziony, Migdałek, Pareek, & Żukowski§, 2020) o en el caso de que no exista un genoma de referencia realizar un ensamble de novo.
Esta carga de trabajo se la realiza con diferentes enfoques de alineación, sin embargo, al igual que cualquier software, cada enfoque presenta diferentes formas de ejecución y consecuentemente de análisis.
Para realizar esta tarea se han desarrollado varios programas como lo son, BWA (Li & Durbin, 2009), Bowtie 2 (JHU, 2016) y otros, cuya función principal es hacer coincidir secuencias cortas provenientes de NGS con secuencias de referencia (Bak, Bodziony, Migdałek, Pareek, & Żukowski§, 2020). Actualmente, BWA es un software muy popular para fines de alineamiento, debido a su compatibilidad con datos de Illumina y la integración de distintos algoritmos de alineamiento con mayor rendimiento para reads de distintas longitudes (Li & Durbin, 2009).
El resultado del mapeo o alineamiento se guarda en un archivo d etexto plano con formato SAM (Sequence Alignment Map) o su equivalente comprimido en binario (BAM). Este archivo contiene de forma codificada toda la información de dónde y cómo cada read se alineo al genoma referencia. De el se pueden luego extraer información de variantes para análisis poblacional o estudios de asociación genómica o información de conteo de reads en regiones codificantes para en análisis de expresión de genes.
El formato SAM consta de un encabezado que comienza con el símbolo @ y una sección de alineamiento que contiene la información de cada uno de los reads que alineo al genoma de referencia. Para conocer más acerca del formato SAM puedes revisar el paper que describe el formato o simplmente en wikipedia. Los archivos SAM se pueden analizar y editar con el software SAMtools. A continuación se muestra un ejemplo del formato SAM con los 11 campos obligatorios del alineamiento y 1 campo opcional.
Formato SAM
BWA
Es un algoritmo de alineación para alinear lecturas de secuencia o secuencias de consulta largas con un genoma de referencia grande como el humano. Elige automáticamente entre alineaciones locales y de extremo a extremo y admite lecturas de extremos emparejados (Li, 2013).
Puedes visualizar las opciones y más información de bwa Aquí
Samtools
Es un conjunto de utilidades que manipulan alineaciones en los formatos SAM (Sequence Alignment Map), BAM y CRAM. Convierte entre los formatos, clasifica, fusiona e indexa, y puede recuperar lecturas en cualquier región rápidamente, puedes tener más información sobre este programa Aquí
Para facilitar el proceso solo se utilizará una muestra la que provieme de la base de datos SRA del NCBI y corresponden a lecturas crudas del salmón del Atlántico Salmo salar en formato fastq, obtenidas por secuenciación de extremos emparejados con un secuenciador Illumina HiSeq2000.
La muestra se alineará al genoma de la mitocondria del salmón del Atlántico puesto que dado que su pequeño tamaño, el alineamiento no debería tardar más de unos pocos minutos.
Descargar Genoma de referencia de la mitocondria e indexar.
Alinear muestra con genoma de referencia.
Analizar el resultado del alineamiento.
Conecte al servidor pomeo usando la siguiente dirección IP: 200.54.220.141 puerto 22.
Para Windows puede conectarse con el software PuTTy y las claves de acceso dadas en la primera clase.
Para configurar el canal bioconda se debe ejecutar el siguiente comando
conda config --add channels bioconda
Instalar los software bwa
conda install -c bioconda bwa
Instalar samtools requiere un proceso de instalación un poco más largo, siga cada paso con atención.
conda install -c bioconda samtools
conda config --add channels bioconda
conda config --add channels conda-forge
conda install samtools==1.11
Si desea constatar la ruta de instalación de un programa determinado puedes colocar el comando “whereis” en la terminal más el programa del que desea obtener la ruta.
whereis sratoolkit
whereis samtools
whereis bwa
La salida de cada comando indicará la ruta de instalación como sigue, para cada programa del que desees obtener información.
/home2/usuario/miniconda3/bin/bwa
Considerando que nuestras secuencias fastq originales tuvieron muy buena calidad, puedes trabajar directamente con ellas. Para ambos casos debes trasladar los archivos a una nueva carpeta, así:
Crear una carpeta denominada en tu usuario de home2
mkdir alineamiento
Ingresa a tu carpeta “alineamiento” y transfiere los archivos de clase anterior colocando los siguientes comandos en la terminal:
mv /home2/usuario/SRA_samples/SRR2006763/SRR2006763_1.fastq /home2/usuario/alineamiento/
mv /home2/usuario/SRA_samples/SRR2006763/SRR2006763_2.fastq /home2/usuario/alineamiento/
Lista tu carpeta de alineamiento para verificar que tienes lo necesario para el alineamiento, hasta ahora deben constar tus dos secuencias “SRR2006763_1.fastq” y “SRR2006763_2.fastq”
ls
Ahora descargaremos el genoma de referencia de la mitocondria de Salmo salar en la misma carpeta de alineamiento, para ello:
En tu navegador entra al siguiente link
Encontrarás una tabla con el listado del genoma de referencia donde se incluyen todos los cromosomas y el genoma de la mitocondria y a cuyo identificardor RefSeq darás clic, como se ve en la siguiente imagen
Listado genoma Salmo salar
Dar clic en la opción FASTA localizada bajo el título e identificador RefSeq de la referencia
Fasta option
Enviar la secuencia FASTA del genoma mitocondrial a un archivo
Genoma mitocondrial
Ve al lugar de descargas de tu computador y encontrarás el archivo denominado “sequence.fasta” al que debes renombrar como “mt.fasta” y subirlo a POMEO, para ello puedes trabajar con WINSCP
Si trabajas con WINSCP debes iniciar sesión, guardar tus datos y conectarte como se indica en la siguiente imagen
WINSCP
Al conectar debes ingresar nuevamente tu contraseña y aparecerá la interfaz de tu servidor con las carpetas que tienes creadas, aquí ingresarás en tu carpeta de alineamiento y arrastrarás el archivo descagado del genoma mitocondrial a la misma.
Alternativamente puedes cargar el genoma iniciando sesión online en el servidor POMEO a traves del puerto 8787 precargado con Rstudio server.
http://200.54.220.141:8787/auth-sign-in
Cuando termines lo anterior puedes ingresar a POMEO, listar tu carpeta de alineamiento con ls y verás tu archivo “mt.fasta” junto con tus secuencias fastq.
Una vez que incluiste en tu carpeta de alineamiento todos los archivos descritos en pasos anteriores podemos proceder a la etapa inicial del alineamiento, que corresponde a la indexación del genoma de referencia con bwa usando el siguiente comando:
bwa index mt.fasta
La salida del comando dará como resultado 5 archivos con extensiones “amb”,“ann”,“bwt”,“pac” y “sa”
Para el alineamiento tendremos las siguientes etapas
Alineamiento de las secuencias contra el genoma de referencia, cuya salida será un archivo con extensión “.sam”
Conversión del archivo SAM a BAM
Inspeccionar el archivo .sam de salida
Ordenar lecturas alineadas por posición
Indexación con Samtools
Exploración de datos con Samtools
Para ejecutar todas las etapas anteriores en ese orden se debe crear un script con nano denominado aln_mt.sh
nano aln_mt.sh
En el script coloca las siguientes instrucciones
#!/bin/bash -l
# para alinear tus dos secuencias fastq al genoma mitocondrial
bwa mem mt.fasta SRR2006763_1.fastq SRR2006763_2.fastq > SRR2006763.sam
# Transformar tu archivo sam a bam
samtools view -Sb -q 30 SRR2006763.sam > SRR2006763.bam
# ordenar tu archivo binario bam
samtools sort SRR2006763.bam -o SRR2006763.sort.bam
# indexar tu archivo bam
samtools index SRR2006763.sort.bam
Ejecuta tu script con bash
bash aln_mt.sh
Ahora que ya tienes tus archivos SAM/BAM puedes observar tu archivo sam con el comando less de linux (recuerda que es un archivo de texto plano)
less SRR2006763.sam
La salida de tu archivo sam se verá como este ejemplo:
Archivo .sam
También puedes realizar análisis adicionales con los siguientes comandos:
| Comando | Descripción |
|---|---|
| samtools flagstat SRR2006763.bam > muestra_stat.txt | análisis estadístico estándar |
| samtools flagstat SRR2006763.bam > muestra_stat.txt | análisis estadístico estándar |
| samtools flags unmap | proporciona los reads no mapeados |
| samtools flags 77 | read 1 - emparejado no mapeado |
| samtools flags 141 | read 2 - emparejado no mapeado |
| samtools flags 99 | Reverse de un read 1 adecuadamente emparejado |
samtools view -f 66 SRR2006763.bam | head -n 10 |
Busca solo reads emparejados en el archivo bam |
Investiga otros comandos de samtools para explorar el alineamiento.
Algunas observaciones adicionales
BWA-MEM se recomienda para secuencias de más de ~70 pb, ya que, generalmente BWA-MEM es más tolerante con los errores dadas las secuencias más largas.
Si quieres ver lo que sucede en cada etapa de alineamiento, puedes colocar en la terminal cada comando del script de forma individual, e.g.
bwa mem mt.fasta SRR2006763_1.fastq SRR2006763_2.fastq > SRR2006763.sam # para alinear tus dos secuencias fastq al genoma mitocondrial samtools view -Sb -q 30 SRR2006763.sam > SRR2006763.bam # Transformar tu archivo sam a bam samtools sort SRR2006763.bam -o SRR2006763.sort.bam # ordenar tu archivo binario bam samtools index SRR2006763.sort.bam # indexar tu archivo bam
El archivo bam no se puede visualizar directamente porque está comprimido y en binario.
Procura listar tu carpeta cada vez que realices un trabajo, una descarga, o ejecutes un script para comprobar tus salidas.
Si quieres comprobar los tamaños de tus archivos al mismo tiempo que los listas, puedes hacerlo con
ls -l -h
Pham-Quoc, C., Kieu-Do, B., & Thinh, T. N. (2019). A high-performance FPGA-based BWA-MEM DNA sequence alignment. Wiley, 1-12.
Li, H., & Durbin, R. (2009). Fast and accurate short read alignment with Burrows-Wheeler Transform. Bioinformatics, 1754-60.
JHU. (2016). Johns Hopkins University. Obtenido de Manual of Bowtie 2 - Fast and sensitive read alignment: http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml